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Abstract A first-order transition is numerically found in a spherical sur- 
face model with skeletons, which are linked to each other at junctions. The 
shape of the triangulated surfaces is maintained by skeletons, which have a 
one-dimensional bending elasticity characterized by the bending rigidity b, 
and the surfaces have no two-dimensional bending elasticity except at the 
junctions. The surfaces swell and become spherical at large b and collapse 
and crumple at small b. These two phases are separated from each other by 
the first-order transition. Although both of the surfaces and the skeleton 
are allowed to self-intersect and, hence, phantom, our results indicate a 
possible phase transition in biological or artificial membranes whose shape 
is maintained by cytoskeletons. 

Keywords Skeleton model ■ Crumpling transition • First-order transition 



1 Introduction 

The so-called hop-diffusion of membrane proteins or lipids observed in cell 
membranes indicates that the cytoskeleton forms compartments on the sur- 
face p] . It has also been recognized that the cytoskeleton of cell membranes 
maintains its cell shape against external forces, and plays important roles 
in cell motion, deformation and in some other functions. The mechanical 
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strength of the cell is considered to be provided by cytoskeletons. Some ar- 
tificial membranes are considered to have skeletons, because they are partly 
polymerized [2]- 

A well-known conventional model for such cell membranes is the cur- 
vature model of Helfrich, Polyakov and Kleinert (HPK) [5P1|5]. The cell 
membranes are considered to be a two-dimensional surface in the HPK pre- 
scription, and therefore the two dimensional differential geometry is suited 
to describe the shape and their corresponding properties [6|[7][8]. However, 
the skeleton is not considered in those curvature models. 

Skeleton models for the cytoskeleton were investigated in [9]. A hard- 
wall and a hard-core potential were assumed on the polymer chains with 
junctions, and the responses to some external stresses and the compression 
modulus were extracted from the Monte Carlo (MC) simulation data [9]. 
The compartmentalized structure was recently investigated in the frame- 
work of the HPK model [TDj . A fluid surface undergoes a phase transition 
due to the compartment. 

However, phase transitions of skeleton surface models are yet to be 
studied. The terminology skeleton surface model in this paper denotes a 
surface model with skeletons; that is surface + skeleton, where the surface 
denotes the ordinary two-dimensional surface and the skeleton denotes one- 
dimensional elastic chains with junctions. It is expected that skeletons can 
make the surface smooth against the surface fluctuations. Therefore it is 
interesting to see whether the crumpling transition occurs in such skeleton 
surface models. The transition is the one that has long been studied theo- 
retically [nHUrrj and numerically P^JT^VT^lfT71[T51fTO|[^gV^T] on the basis 
of the HPK model, and an experimental investigation on the transition has 
also been performed recently [2]. 

We consider that the problem of whether the skeleton surface model 
undergoes the phase transition or not is non-trivial. In fact, we know that 
one-dimensional objects, such as an elastic ring obeying a Hamiltonian for 
the local interactions, has no phase transition. 

In this paper, we study the phase structure of a spherical skeleton sur- 
face model, which is defined on a triangulated spherical surface with poly- 
mer chains linked to each other at junctions, focusing on how the smooth 
phase separated from the crumpled phase. The Hamiltonian is a linear com- 
bination of the Gaussian bond potential Si, the one-dimensional bending 
energy 5*2 on the polymer chain, and the two-dimensional bending energy 
Sj at the junctions. 

The interaction in the model is not purely one-dimensional. One-dimensional 
polymer chains have two-dimensional interactions through the two-dimensional 
potential S\ and the two-dimensional bending energy Sj at the junctions, 
which are considered to be the boundary of polymer chains. However, it 
is obvious that the elastic skeletons maintain the shape of the membrane, 
because triangulated surfaces have no two-dimensional bending elasticity 
except at the junctions of skeletons. 
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We will find a first-order transition in the model under certain con- 
ditions on the length of the polymer chains and the bending rigidity for 
Sj. Although the model is allowed to self-intersect and, hence, phantom, 
our results indicate a possible phase transition in biological or artificial 
membranes whose mechanical strength is provided by the cytoskeletons. 



2 Model 

The triangulated surfaces are characterized by N the total number of ver- 
tices, Ns the total number of vertices on the chains, Nj the total num- 
ber of junctions, and L the length of chains between junctions. Figures 
[l^a) andQIb) show surfaces of size (N,N s ,Nj,L) = (2562,600,42,6) and 
(N, N s , Nj, L) = (4002, 840, 42, 8), respectively. They are assumed to be the 
starting configurations of the MC simulations. Thick lines denote the com- 
partment boundary, which forms chains and the junctions. All the vertices 
can fluctuate only locally, and they are prohibited from diffusing because 
of the fixed connectivity nature of the lattice. 




(a) (2562, 600, 42, 6) (b) (4002, 840, 42, 8) 



Fig. 1 Starting configurations of (a) [N, Ns, Nj,L) = (2562, 600, 42, 6) and (b) 
(TV, N s , Nj, L) = (4002, 840, 42, 8), where N is the total number of vertices, N s 
is the total number of vertices on the chains, Nj is the total number of junc- 
tions, and L is the length of chains between junctions. Thick lines denote the 
compartment boundary, which forms chains and junctions. 

A triangulated surface of size N — 10£ 2 + 2 (= the total number of 
vertices) is obtained by dividing every edge of the icosahedron into £ pieces 
of uniform length; I is the edge length of triangles in the icosahedron. The 
configurations are thus characterized by N 5 = 12 and iV 6 = N~ 12, where N q 
is the total number of vertices with co-ordination number q. Compartment 
structures are obtained by dividing £ further into m pieces (m = 1, 2, • • •), 
and we have chains of the uniform length L=(£/m)—2. The reason for the 
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subtraction of —2 is because of the junctions at the two end points of the 
chain. On the surfaces in Figs. UX&) and [lib), we have L = 6, I — 16 and 
L = 8, £ — 20, respectively. 

The compartment structures shown in Figs. HJa) and HJb) are identi- 
cal to those from the model in Ref. [10] . where the boundary bonds of 
the compartment were kept unflipped in the dynamical triangulation. The 
boundary of compartment in Ref. jTD] is identical with the linear chains with 
junctions in this paper. The junctions form hexagons or pentagons; the to- 
tal number of pentagon is 12 and the remaining junctions are hexagons. 
The total number of the compartments depends on the surface size N, and 
therefore, it is increased with increasing TV. However, the chain length L 
can be chosen to be constant and independent of N. We fix the chain length 
L to the following two values: 

L = 6, L = 8, (1) 

which respectively correspond to the values n = 21, n = 36, the total number 
of vertices inside a compartment [lOj . The reason why we fix n is that the 
size of compartment is considered to be finite and, thus, it is expected 
that total number of lipids in the compartment also remains finite in the 
cell membranes. Note also that n does not exactly correspond to the total 
number of lipids in a compartment of cell membranes. 

The Hamiltonian of the model is given by a linear combination of the 
Gaussian bond potential Si, the one-dimensional bending energy £2, and 
the two-dimensional bending energy Sj, which are defined by 

5i = E - x rf < 52 = E(! - cos %))> 

S J = 5^(1 -n< ■!!,•), (2) 

m 

where Y^Uj) m $1 is the sum over the bond (ij) connecting the vertices i 
and j, and m ^2 is the sum over bonds i and j, which contain not only 

bonds in the chains but also bonds that connect the center of the junction 
and the neighboring vertices of the chains. Si is defined not only on the 
skeleton but also on the surface, while S2 is defined only on the chains. 
y)/ij\ in Sj is the sum over triangles i and j, which share the central point 
of the junction as the common vertex and form a hexagonal or pentagonal 
junction. The symbol Bm) in S2 is the angle between the bonds i and j, 
and rii in Sj is the unit normal vector of the triangle i at the junctions. 

Figure [2] shows a hexagonal junction linked to chains, the unit normal 
vectors n,;, rij, and the angle 0(ij), which is defined not only on the vertices 
of chains but also on those nearest to the center of the junction. 

We note that Si includes bonds which are the boundary edges of the 
hexagons (or the pentagons). The reason of this inclusion is for the sake 



Fig. 2 A hexagonal junction linked to chains. The unit normal vector n; in Sj 
is defined on the triangle i in the hexagon, and the angle in 5*2 is defined not 
only on the vertices of chains but also on those nearest to the center of junction. 

of an in-plane elasticity of surfaces at the junctions. If these bonds were 
excluded from Si, chains would be able to freely move into the in-plane 
directions at the junctions. 

The partition function Z of the model is defined by 

( l N 

Z= Y[dX ie xp[-S{X)}, (3) 

»=i 

S(X) = S 1 + bS 2 + bjS J , 

where b is the one-dimensional bending rigidity, bj is the two-dimensional 
bending rigidity at the junctions, and J denotes that the center of the 
surface is fixed. In this paper, bj is fixed to 

bj = 10, (4) 

so that the junctions are sufficiently smooth. The value bj = 10 is relatively 
larger than the first-order transition point 6 C ~0.8 in the tethered surface 
model [T6lll7j . Therefore, the hexagonal or pentagonal junctions are almost 
flat even when the surface is in the crumpled phase at sufficiently small b. 

Both b and bj have units of kT, where k is the Boltzmann constant, and 
T is the temperature. The surface tension coefficient a of Si is fixed to a = 1; 
this is always possible because of the scale invariant property of the model. 
In fact, from the expression aSi+bS2+bjSj we immediately understand that 
a = 1 is possible, because the factor a of Si can be eliminated due to the 
scale invariance of the partition function. Since the unit of a is (1/length) 2 , 
the length unit of the model is given by \flja. We use the unit of length 
provided by yjl/a=l in this paper, although a can arbitrarily be chosen. 



6 



The difference between the model in this paper and the conventional 
tethered surface model such as the one in [16] must be emphasized. The sur- 
face shape is maintained by the two-dimensional bending energy defined all 
over the surface in the conventional surface model. On the contrary, the sur- 
face shape of the model in this paper is maintained by the one-dimensional 
S2 and two-dimensional Sj; both are defined only on the skeletons. 



3 Monte Carlo technique 

The canonical Metropolis technique is used to update the variable X. The 
update of X is divided into two steps. The first is the update of X in the 
chains: the X are shifted so that X' = X+SX, where SX is randomly chosen 
in a small sphere. The new position X' is accepted with the probability 
Min[l, exp(— AS)), where AS = 5(new) — 5 (old). The second step involves 
the update of X in the junctions, and this is further divided into three 
processes: the first is a random and simultaneous shift of 7 (or 6) vertices of 
the hexagon (or the pentagon) including the central vertex, and the second 
and the third are a random translation and a random rotation of these 
vertices by assuming the hexagon (or the pentagon) to be a rigid object. 
All of these MC processes are performed under about a 50% acceptance 
rate, which is controlled by small numbers fixed at the beginning of the 
simulations. We introduce the lower bound of 1 x 1CP 8 for the area of 
triangles. No lower bound is imposed on the bond length. 

We use surfaces of size (N, N s , Nj, L) = (2562, 600, 42, 6), (5762, 1350, 92, 6), 
(10242,2400,162,6), and (16002,3750,252,6) for the length L = 6, and 
(N, N s , Nj, L) = (4002, 840, 42, 8), (9002, 1890, 92, 8), and (16002, 3360, 162, 8) 
for the length L = 8. A random number sequence called Mersenne Twister 
l22l is used in the simulations. 



4 Results 

First, we show snapshots of surface obtained at b — 8.1 (smooth phase) 
and 6 = 8 (crumpled phase) in FigsE^a) and[3jb), respectively. The size 
of surface is (N, N s , Nj, L) = (16002, 3750, 252, 6). The surface sections are 
shown in Figs [3jc) l3^d) . We immediately find that the surface is smooth 
(crumpled) at 6 = 8.1 (6 = 8). 

The surface size is reflected in the mean square size of X 2 , which is 
defined by 

i i 

where X is the center of the surface. We have X 2 ~ 101 and X 2 ~ 10.1 on 
the surfaces in Figs. [3fa) and 0Kb), respectively. 
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(a) Smooth surface at 6 = 8.1 (b) Collapsed surface at 6 = 8 




(c) The surface section (d) The surface section 



Fig. 3 Snapshot of the surface of size (N, N s ,Nj, L) = (16002, 3750, 252, 6) ob- 
tained in the smooth phase at (a) 6 = 8.1 and in the crumpled phase at (b) 6 = 8, 
both of which are close to the transition point. The mean square sizes of X 2 are 
(a) X 2 ~ 101 and (b) X 2 ~ 10.1. 



The Gaussian bond potential Si/N against b is shown in Figs. [H a) and 
0Jb), where the length L of chains between junctions are L = 6 and L = 8, 
respectively. The surface size is N = 5762, N — 10242, and N = 16002 for 
L = 6, and N = 4002, N = 9002, and N = 16002 for L = 8. The solid lines 
connecting the data are drawn by a multihistogram reweighting technique 
[23j . Because of the scale invariant property of the partition function, we 
have Si/N = 3(N - 1)/2N ~ 1.5. All of the results presented in Figs. 
SJ a )HKt>) satisfy the expected relation. 

The mean square size of X 2 versus b is plotted in Figs. OJa) and[UJb), 
where L = 6 and L = 8, respectively. The data was obtained on the same 
sized surfaces as those in Figs.|4}a) and[4jb). The solid lines drawn on the 
data were obtained by a multihistogram reweighting technique. We find 
a sharp change of X 2 both in Figs. EJa) and EJb) as N increases. This 
sharp change of X 2 seems discontinuous and hence indicates a first-order 
transition. We find also that the transition point b c moves left on the 6-axis 
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Fig. 4 The Gaussian bond potential Si/N vs. b obtained on the surfaces of (a) 
L = 6 and (b) L — 8. Si/N is predicted to be Si/N ~ 1.5. The two-dimensional 
bending rigidity bj at the junctions was fixed to bj = 10. 



as N increases in both cases (i.e., L — 6 and L — 8). This implies that the 
finite-size effect still remained on large surfaces of size N = 16002. However, 
we should recall that the bending energy S2 is defined on the chains of size 
Ns = 3750 and N$ = 3360, which are relatively small, on such large surfaces. 
It is expected that more and more large surfaces are necessary to remove 
the finite-size effect. We find also that the surface softens as the length of 
chain increases. In fact, the transition point b c moves right on the &-axis as 
L increases from L = 6 to L = 8, as can be seen in Figs.[5ja) and[5^b). 




Fig. 5 The mean square size of X 2 versus b is obtained on the surfaces of (a) 
L = 6 and (b) L = 8. The rigidity bj was fixed to bj = 10. The curves are drawn 
by the multihistogram reweighting technique. 
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It must be noted that a large surface in the collapsed phase hardly 
turns into the smooth phase even at the transition point. This phenomenon 
seems typical to surface simulations [T6lfT7] based on the canonical MC 
simulation technique. When a large surface configuration is once trapped 
in a potential minimum, the configuration appears almost confined inside 
the potential valley. These troublesome phenomena can be avoided with 
more sophisticated MC techniques [241125] , 

We also comment on the Hausdorff dimension H , which is defined by the 
scaling relation X 2 oc N 2 / H and was not calculated explicitly. However, it 
is obvious that H ~ 2 is satisfied in the smooth phase close to the transition 
point, because the spherical surfaces are swollen and smooth at that point 
as can be seen in the snapshot of Figs. [3Ja) and[3]Jc). Moreover, we expect 
that H is greater than the physical bound H > 3 in the crumpled phase 
close to the transition point. In fact, we find from Figs. 02a) and E^b) 
that X 2 appears to remain unchanged as N increases in the crumpled 
phase. This implies that H is very large and is in contrast to that of the 
conventional tethered surface model in p~6|.l 1 T] . where H is less than the 
physical bound in the crumpled phase. 

The bending energy obtained on the surfaces of L = 6 and L = 8 

is shown in Figs. HJa) and [S^b), where N' s — N$ + 6Nj — 12 is the total 
number of vertices on which S2 is defined. We have N' s = 5250 for the 
surface of (N, N s ,Nj, L) = (16002, 3750, 252, 6) and N' s = 4320 for that of 
(N, N s , Nj, L) = (16002, 3360, 162, 8). Although the data in the figures are 
connected with smooth curves obtained by a multihistogram reweighting 
technique, the discontinuous nature of £2/-/^ is almost apparent. 




Fig. 6 The one-dimensional bending energy Ss/N'g vs. b obtained on the surfaces 
of (a) L = 6 and (b) L = 8. 
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The phase transition is not reflected in the two-dimensional bending 
energy Sj. Since the corresponding bending rigidity bj was chosen to make 
bj = 10 (as shown in Eq.((4])), the junctions are maintained smooth enough. 
Therefore, Sj defined on the junctions is kept small and exhibits only weak 
fluctuations. 

On the other hand, it is possible to define a two-dimensional bending 
energy 

4 2) =^(l-n 4 .n,) (6) 

(«» 

on the surface, where is the unit normal vector of the triangle i, and 
J2{{ij)) i s tne sum over au nearest neighbor triangles {(ij)) on the surface 

(2) 

except at the junctions. Although S 2 is not included in the Hamiltonian, 

the transition can be reflected in . Figures[7{a) and[7Ub) show S 2 /Nb 
against b obtained on the surfaces of L = 6 and L — 8, where Nb is the total 

(2) (2) 

number of bonds where S 2 is defined. Note that the bonds where S 2 is 
undefined are those where Sj is defined. We immediately find the expected 
behavior of S 2 , which again indicates a discontinuous transition. 



S 2 lz 7N B S 2 lz 7N : 




7.8 8 , 8.2 8.4 13 , 14 

(a) b (b) b 

Fig. 7 The two-dimensional bending energy \ /Nb against b obtained on the 
surfaces of (a) L = 6 and (b) L = 8. is defined by Eq.© and is not included 
in the Hamiltonian. Nb is the total number of bonds where S% is defined. 



The specific heat for S2 defined by 



Cs 2 = -£r((S2~(S 2 )) 2 ) (7) 
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has an anomalous behavior if the model has the phase transition. Figures 
[HJa) and EJb) show Cs 2 against b obtained on the surfaces of L = 6 and 
L = 8. Sharp peaks of Cs 2 are considered to be a sign of the phase transition. 




Fig. 8 The specific heat Cs 2 for S2 against b obtained on the surfaces of (a) 
L — 6 and (b) L — S. Cs 2 is defined by Eq.([7J). The error bars are the statistical 
error, which is obtained by the binning analysis. 

(2) 

Another specific heat, which is the variance of S% , can also be defined 
such that 

C s ™=jjj((S™-(S™)) 2 ), (8) 

where N'=N— 6Nj+12 is the total number of vertices N minus 6Nj — 12 
the number of vertices nearest to the junctions. Figures [9]^a) and[9^b) show 
C„(2) against b obtained on the surfaces of L = 6 and L = 8. We find that 

C„(2) has the same anomalous behavior as that of Cs 2 in Figs. EJa) and 

Mb)- 

In order to see the scaling property of Cs 2 and C„(2) , we plot the peak 

values Cg ax against N' s and C m ( tf against N' in Figs. HD^a) and[lOi;b), 

respectively, in a log-log scale. The error bars on the peak values are the 
statistical errors, which were obtained by the so-called binning analysis. 
The straight lines were drawn by fitting the three largest data values to 

cr«(4f, c;^(x(7vr 2 , (9) 

where 01, 02 are critical exponents. Thus, we have 

o-i = 1.27 ±0.19, 0-2 = 1.34 ±0.11, (X = 6), 

01 = 1.16 ±0.13, 0-2 = 1.12 ± 0.05, (L = 8). (10) 
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Fig. 9 The specific heat C (2) for So 2 ' against b obtained on the surfaces of (a) 
L = 6 and (b) L = 8. C (2) is defined by Eq.®. The error bars are the statistical 
error, which is obtained also by the binning analysis. 



These values are slightly larger than a = 1, however, they are consistent 
with the discontinuous transition, which was already expected by all of the 
above-presented results. 



5 Summary and Conclusion 

We have investigated phase transitions between the smooth phase and the 

crumpled phase of a surface model with elastic chains joined to each other 

at junctions by using the canonical Monte Carlo simulation technique. The 

model is defined on uniformly triangulated spheres, which were obtained by 

splitting the triangles of the icosahedron. A compartmentalized structure 

was built on the surface, and the boundary of the compartment forms 

the skeletons, which have a one-dimensional bending energy and maintain 

shape of the surface. No two-dimensional bending energy is assumed on 

the surface except at the junctions. We assumed a high elasticity at the 

junctions; bj — 10, and that the length L of chains between junctions is 

relatively small; L = 6 and L = 8. 

A hrst-order transition has been found. The bending energy £2 of the 

(2) 

chains and the two-dimensional bending energy , which is not included 
in the Hamiltonian, are found to have a jump at finite b. The corresponding 
specific heats have an anomalous peak, which is typical of a first-order 
transition. The mean square size X 2 also has a clear jump at the transition 
point, indicating that the transition distinguishes the swollen (or smooth) 
phase from the collapsed (or crumpled) phase. 
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on the surfaces of L = 6 and L = 8. The straight lines are drawn by fitting the 
largest three data of C™^ and C m ( 2* to Eq. © . The peak values and the statistical 

errors were obtained by multihistogram reweighting. 



The fact that a skeleton surface model has a phase transition is re- 
markable because it is well-known that we can see no phase transition in 
one-dimensional objects which are governed by a Hamiltonian for the local 
interactions. The elasticity of the junctions are considered to be a reason 
for the phase transition. The chains in our model share the two-dimensional 
Gaussian bond potential Si and the two-dimensional bending energy Sj at 
the boundaries, which are the junctions. However, the two-dimensional Si 
seems to play no significant role in the transition, because a skeleton sur- 
face model with a one-dimensional bond potential also undergoes a phase 
transition [55]. It is expected that the skeleton surface model has a phase 
transition even without Si as in HPK model [IT] . 

The phase structure of the skeleton surface model seems to be signif- 
icantly dependent on the junctions. Therefore, a skeleton surface model 
with rigid junctions is interesting. A rigid junction is a junction that has 
not only the infinite bending elasticity but also the infinite in-plane elastic- 
ity and, consequently, it is different from the junctions in this paper, which 
have finite in-plane elasticity even in the limit of bj— >oo. The fluid nature 
can also be incorporated in the skeleton surface model with the dynamical 
triangulation MC technique, and the phase structure of the model seems 
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to be changed. Many interesting points remain to be studied in skeleton 
surface models. 
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